Data and code for paper on evolution during a 2017 epidemic of Pasteuria ramosa in Little Appleton Lake. Authors: Camden D. Gowler Haley Essington Patrick A. Clay Bruce O’Brien Clara L. Shaw Rebecca Bilich Meghan A. Duffy

Analysis of data collected on Little Appleton’s 2017 Pasteuria epidemic

# Calculating the proportion that never reproduced
LA_infected <- LA_infected %>%
  mutate(anyrepro = if_else(clutches == "0", "0", "1"))

LA_infected$anyrepro <- as.numeric(LA_infected$anyrepro)

# Adding a column for log spores
LA_infected$logspores <- log10(LA_infected$spores_tot)
LA_infected$lnspores <- log(LA_infected$spores_tot)
LA_infected <- subset(LA_infected, host_time_point != "4")
str(LA_infected)
## tibble[,16] [259 × 16] (S3: tbl_df/tbl/data.frame)
##  $ clone          : num [1:259] 343 17 304 304 324 311 323 29 333 29 ...
##  $ block          : num [1:259] 1 2 1 1 2 1 1 1 2 1 ...
##  $ rep            : chr [1:259] "R01" "R2" "R03" "R02" ...
##  $ treatment      : chr [1:259] "B" "C" "B" "C" ...
##  $ death_date     : Date[1:259], format: "2019-02-18" "2019-02-22" ...
##  $ infection      : chr [1:259] NA NA NA NA ...
##  $ birth_date     : Date[1:259], format: "2019-02-09" "2019-02-14" ...
##  $ host_time_point: num [1:259] 3 1 3 3 3 3 3 1 3 1 ...
##  $ para_time_point: chr [1:259] "3" "none" "3" "none" ...
##  $ lifespan       : num [1:259] 9 8 15 15 12 22 22 22 17 23 ...
##  $ avg            : num [1:259] 6.75 0 0.25 0 0 ...
##  $ spores_tot     : num [1:259] 6750 0 250 0 0 ...
##  $ clutches       : num [1:259] 0 0 4 2 0 0 0 2 0 1 ...
##  $ anyrepro       : num [1:259] 0 0 1 1 0 0 0 1 0 1 ...
##  $ logspores      : num [1:259] 3.83 -Inf 2.4 -Inf -Inf ...
##  $ lnspores       : num [1:259] 8.82 -Inf 5.52 -Inf -Inf ...
LA_infected$clone <- as.factor(as.character(LA_infected$clone))

# get mean lifespan for each clone*treatment
stats_lifesp <- LA_infected %>%
  group_by(host_time_point, para_time_point, clone) %>%
  summarise(mn_lifespan = mean(lifespan)) 
## `summarise()` has grouped output by 'host_time_point', 'para_time_point'. You can override using the `.groups` argument.
# for selecting contemporaneous pairs of hosts/parasites
stats_lifesp$host_time_point <- as.character(stats_lifesp$host_time_point)
stats_lifesp$para_time_point <- as.character(stats_lifesp$para_time_point)

LA_inf <- stats_lifesp %>% # select non-controls
  filter(para_time_point != "none") 

LA_control <- stats_lifesp %>% # select controls
  filter(host_time_point != "4") %>%
  filter(para_time_point == "none") %>%
  ungroup() %>%
  select(-para_time_point)

colnames(LA_control)[3] <- "control_lifespan"

# calculate average lifespan for all controls together (grouped by time point)
LA_control_avg <- LA_control %>%
  group_by(host_time_point) %>%
  summarise(mean_life = mean(control_lifespan))

# turn into an unamed matrix so we can easily extract values later
mean_mat <- unname(as.matrix(LA_control_avg))

# rejoin infected and control dfs
LA_rel_vir_lifespan <- full_join(LA_inf, LA_control)
## Joining, by = c("host_time_point", "clone")
LA_rel_vir_lifespan$host_time_point <- as.factor(LA_rel_vir_lifespan$host_time_point)

# make a separate file to use later for infected v. control comparison
LA_rel_vir_lifespan_infvcontrol <- full_join(LA_inf, LA_control)
## Joining, by = c("host_time_point", "clone")
LA_rel_vir_lifespan_infvcontrol$host_time_point <- as.factor(LA_rel_vir_lifespan_infvcontrol$host_time_point)

# Replace NAs with the average control lifespan for that group
# do this b/c didn't have controls for each clone
LA_rel_vir_lifespan$control_lifespan[is.na(LA_rel_vir_lifespan$control_lifespan) & 
                              LA_rel_vir_lifespan$host_time_point == "1"] <- mean_mat[1,2] # means
LA_rel_vir_lifespan$control_lifespan[is.na(LA_rel_vir_lifespan$control_lifespan) & 
                              LA_rel_vir_lifespan$host_time_point == "2"] <- mean_mat[2,2]
LA_rel_vir_lifespan$control_lifespan[is.na(LA_rel_vir_lifespan$control_lifespan) & 
                              LA_rel_vir_lifespan$host_time_point == "3"] <- mean_mat[3,2]

# calculate relative virulence using mean and median control lifespans
LA_rel_vir_lifespan$control_lifespan <- as.numeric(LA_rel_vir_lifespan$control_lifespan)

LA_rel_vir_lifespan <- LA_rel_vir_lifespan %>%
  mutate(rel_lifespan_mn = mn_lifespan/control_lifespan)

LA_rel_vir_lifespan <- as.data.frame(LA_rel_vir_lifespan)

str(LA_rel_vir_lifespan)
## 'data.frame':    41 obs. of  6 variables:
##  $ host_time_point : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 2 2 ...
##  $ para_time_point : chr  "1" "1" "1" "1" ...
##  $ clone           : Factor w/ 28 levels "17","211","22",..: 1 3 13 16 2 4 5 6 7 8 ...
##  $ mn_lifespan     : num  33 38.7 27.1 42 35.8 ...
##  $ control_lifespan: num  46 37.7 58.8 26 54.9 ...
##  $ rel_lifespan_mn : num  0.717 1.027 0.461 1.615 0.651 ...
str(LA_infected)
## tibble[,16] [259 × 16] (S3: tbl_df/tbl/data.frame)
##  $ clone          : Factor w/ 28 levels "17","211","22",..: 28 1 14 14 21 17 20 13 24 13 ...
##  $ block          : num [1:259] 1 2 1 1 2 1 1 1 2 1 ...
##  $ rep            : chr [1:259] "R01" "R2" "R03" "R02" ...
##  $ treatment      : chr [1:259] "B" "C" "B" "C" ...
##  $ death_date     : Date[1:259], format: "2019-02-18" "2019-02-22" ...
##  $ infection      : chr [1:259] NA NA NA NA ...
##  $ birth_date     : Date[1:259], format: "2019-02-09" "2019-02-14" ...
##  $ host_time_point: num [1:259] 3 1 3 3 3 3 3 1 3 1 ...
##  $ para_time_point: chr [1:259] "3" "none" "3" "none" ...
##  $ lifespan       : num [1:259] 9 8 15 15 12 22 22 22 17 23 ...
##  $ avg            : num [1:259] 6.75 0 0.25 0 0 ...
##  $ spores_tot     : num [1:259] 6750 0 250 0 0 ...
##  $ clutches       : num [1:259] 0 0 4 2 0 0 0 2 0 1 ...
##  $ anyrepro       : num [1:259] 0 0 1 1 0 0 0 1 0 1 ...
##  $ logspores      : num [1:259] 3.83 -Inf 2.4 -Inf -Inf ...
##  $ lnspores       : num [1:259] 8.82 -Inf 5.52 -Inf -Inf ...
# get mean number clutches for each clone*treatment
stats_clutches <- LA_infected %>%
  group_by(host_time_point, para_time_point, clone) %>%
  summarise(mn_clutches = mean(clutches)) 
## `summarise()` has grouped output by 'host_time_point', 'para_time_point'. You can override using the `.groups` argument.
# for selecting contemporaneous pairs of hosts/parasites
stats_clutches$host_time_point <- as.character(stats_clutches$host_time_point)
stats_clutches$para_time_point <- as.character(stats_clutches$para_time_point)

LA_inf <- stats_clutches %>% # select non-controls
  filter(para_time_point != "none") 

LA_control <- stats_clutches %>% # select controls
  filter(host_time_point != "4") %>%
  filter(para_time_point == "none") %>%
  ungroup() %>%
  select(-para_time_point)

colnames(LA_control)[3] <- "control_clutches"

# calculate average no. clutches for all controls together (grouped by time point)
LA_control_avg <- LA_control %>%
  group_by(host_time_point) %>%
  summarise(mean_clutches = mean(control_clutches))

# turn into an unamed matrix so we can easily extract values later
mean_mat <- unname(as.matrix(LA_control_avg))

# rejoin infected and control dfs
LA_rel_vir_clutches <- full_join(LA_inf, LA_control)
## Joining, by = c("host_time_point", "clone")
LA_rel_vir_clutches$host_time_point <- as.factor(LA_rel_vir_clutches$host_time_point)

# make a separate file to use later for infected v. control comparison
LA_rel_vir_clutches_infvcontrol <- full_join(LA_inf, LA_control)
## Joining, by = c("host_time_point", "clone")
LA_rel_vir_clutches_infvcontrol$host_time_point <- as.factor(LA_rel_vir_clutches_infvcontrol$host_time_point)

# Replace NAs with the average control number of clutches for that group
# do this b/c didn't have controls for each clone
LA_rel_vir_clutches$control_clutches[is.na(LA_rel_vir_clutches$control_clutches) & 
                              LA_rel_vir_clutches$host_time_point == "1"] <- mean_mat[1,2] # means
LA_rel_vir_clutches$control_clutches[is.na(LA_rel_vir_clutches$control_clutches) & 
                              LA_rel_vir_clutches$host_time_point == "2"] <- mean_mat[2,2]
LA_rel_vir_clutches$control_clutches[is.na(LA_rel_vir_clutches$control_clutches) & 
                              LA_rel_vir_clutches$host_time_point == "3"] <- mean_mat[3,2]

# calculate relative virulence using mean and median control clutches
LA_rel_vir_clutches$control_clutches <- as.numeric(LA_rel_vir_clutches$control_clutches)

LA_rel_vir_clutches <- LA_rel_vir_clutches %>%
  mutate(rel_clutches_mn = mn_clutches/control_clutches)

LA_rel_vir_clutches <- as.data.frame(LA_rel_vir_clutches)

str(LA_rel_vir_clutches)
## 'data.frame':    41 obs. of  6 variables:
##  $ host_time_point : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 2 2 ...
##  $ para_time_point : chr  "1" "1" "1" "1" ...
##  $ clone           : Factor w/ 28 levels "17","211","22",..: 1 3 13 16 2 4 5 6 7 8 ...
##  $ mn_clutches     : num  0 5.333 0.867 0 0.25 ...
##  $ control_clutches: num  11.3 9 19.4 1 16 ...
##  $ rel_clutches_mn : num  0 0.5926 0.0447 0 0.0156 ...
LA_rel_vir_2 <- full_join(LA_rel_vir_lifespan, LA_rel_vir_clutches)
## Joining, by = c("host_time_point", "para_time_point", "clone")
str(LA_rel_vir_2)
## 'data.frame':    41 obs. of  9 variables:
##  $ host_time_point : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 2 2 ...
##  $ para_time_point : chr  "1" "1" "1" "1" ...
##  $ clone           : Factor w/ 28 levels "17","211","22",..: 1 3 13 16 2 4 5 6 7 8 ...
##  $ mn_lifespan     : num  33 38.7 27.1 42 35.8 ...
##  $ control_lifespan: num  46 37.7 58.8 26 54.9 ...
##  $ rel_lifespan_mn : num  0.717 1.027 0.461 1.615 0.651 ...
##  $ mn_clutches     : num  0 5.333 0.867 0 0.25 ...
##  $ control_clutches: num  11.3 9 19.4 1 16 ...
##  $ rel_clutches_mn : num  0 0.5926 0.0447 0 0.0156 ...
str(LA_infected)
## tibble[,16] [259 × 16] (S3: tbl_df/tbl/data.frame)
##  $ clone          : Factor w/ 28 levels "17","211","22",..: 28 1 14 14 21 17 20 13 24 13 ...
##  $ block          : num [1:259] 1 2 1 1 2 1 1 1 2 1 ...
##  $ rep            : chr [1:259] "R01" "R2" "R03" "R02" ...
##  $ treatment      : chr [1:259] "B" "C" "B" "C" ...
##  $ death_date     : Date[1:259], format: "2019-02-18" "2019-02-22" ...
##  $ infection      : chr [1:259] NA NA NA NA ...
##  $ birth_date     : Date[1:259], format: "2019-02-09" "2019-02-14" ...
##  $ host_time_point: num [1:259] 3 1 3 3 3 3 3 1 3 1 ...
##  $ para_time_point: chr [1:259] "3" "none" "3" "none" ...
##  $ lifespan       : num [1:259] 9 8 15 15 12 22 22 22 17 23 ...
##  $ avg            : num [1:259] 6.75 0 0.25 0 0 ...
##  $ spores_tot     : num [1:259] 6750 0 250 0 0 ...
##  $ clutches       : num [1:259] 0 0 4 2 0 0 0 2 0 1 ...
##  $ anyrepro       : num [1:259] 0 0 1 1 0 0 0 1 0 1 ...
##  $ logspores      : num [1:259] 3.83 -Inf 2.4 -Inf -Inf ...
##  $ lnspores       : num [1:259] 8.82 -Inf 5.52 -Inf -Inf ...
# get mean number clutches for each clone*treatment
stats_spores <- LA_infected %>%
  group_by(host_time_point, para_time_point, clone) %>%
  summarise(mn_spores = mean(spores_tot)) 
## `summarise()` has grouped output by 'host_time_point', 'para_time_point'. You can override using the `.groups` argument.
# for selecting contemporaneous pairs of hosts/parasites
stats_spores$host_time_point <- as.character(stats_spores$host_time_point)
stats_spores$para_time_point <- as.character(stats_spores$para_time_point)

LA_inf <- stats_spores %>% # select non-controls
  filter(para_time_point != "none") 

# join spore data with earlier calculations
LA_means <- full_join(LA_rel_vir_2, LA_inf)
## Joining, by = c("host_time_point", "para_time_point", "clone")
LA_means$host_time_point <- as.factor(LA_means$host_time_point)

str(LA_means)
## 'data.frame':    41 obs. of  10 variables:
##  $ host_time_point : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 2 2 ...
##  $ para_time_point : chr  "1" "1" "1" "1" ...
##  $ clone           : Factor w/ 28 levels "17","211","22",..: 1 3 13 16 2 4 5 6 7 8 ...
##  $ mn_lifespan     : num  33 38.7 27.1 42 35.8 ...
##  $ control_lifespan: num  46 37.7 58.8 26 54.9 ...
##  $ rel_lifespan_mn : num  0.717 1.027 0.461 1.615 0.651 ...
##  $ mn_clutches     : num  0 5.333 0.867 0 0.25 ...
##  $ control_clutches: num  11.3 9 19.4 1 16 ...
##  $ rel_clutches_mn : num  0 0.5926 0.0447 0 0.0156 ...
##  $ mn_spores       : num  41000 46750 166750 263750 260625 ...

Infected vs. uninfected comparison

Let’s start by looking at lifespan of control/unexposed animals vs. infected animals. Based on what we know about the system, these may not be very different.

# Making a long data set with the uninfected & infected animals to compare lifespan
LA_rel_vir_lifespan_infvcontrol <- LA_rel_vir_lifespan_infvcontrol %>%
  rename(Infected = mn_lifespan)

LA_rel_vir_lifespan_infvcontrol <- LA_rel_vir_lifespan_infvcontrol %>%
  rename(Unexposed = control_lifespan)

LA_rel_vir_lifespan_infvcontrol_long <- LA_rel_vir_lifespan_infvcontrol %>%
  gather(infvcontrol, lifespan, Infected:Unexposed)

LA_rel_vir_lifespan_infvcontrol_long <- na.omit(LA_rel_vir_lifespan_infvcontrol_long)
  

# Now do the same for the reproduction data
LA_rel_vir_clutches_infvcontrol <- LA_rel_vir_clutches_infvcontrol %>%
  rename(Infected = mn_clutches)

LA_rel_vir_clutches_infvcontrol <- LA_rel_vir_clutches_infvcontrol %>%
  rename(Unexposed = control_clutches)

LA_rel_vir_clutches_infvcontrol_long <- LA_rel_vir_clutches_infvcontrol %>%
  gather(infvcontrol, clutches, Infected:Unexposed)

LA_rel_vir_clutches_infvcontrol_long <- na.omit(LA_rel_vir_clutches_infvcontrol_long)
lifespaninfvcontrol <- ggplot(LA_rel_vir_lifespan_infvcontrol_long, 
                               aes(x=infvcontrol, y=lifespan)) + 
  geom_violin() + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33') +
  xlab("Infected or unexposed") + 
  ylab("Mean lifespan \n(days)")  +
  theme_cowplot()

lifespaninfvcontrol

infvcontrollifespanmodel <- glm(lifespan ~ infvcontrol, 
             family = "quasipoisson", 
             data = LA_rel_vir_lifespan_infvcontrol_long)

plot(infvcontrollifespanmodel) 

summary(infvcontrollifespanmodel)
## 
## Call:
## glm(formula = lifespan ~ infvcontrol, family = "quasipoisson", 
##     data = LA_rel_vir_lifespan_infvcontrol_long)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2551  -0.3763  -0.0037   0.5306   1.7622  
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.55025    0.02475 143.435  < 2e-16 ***
## infvcontrolUnexposed  0.28494    0.04235   6.729 1.05e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.8746756)
## 
##     Null deviance: 87.448  on 56  degrees of freedom
## Residual deviance: 49.126  on 55  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
tapply(LA_rel_vir_lifespan_infvcontrol_long$lifespan, LA_rel_vir_lifespan_infvcontrol_long$infvcontrol, mean)
##  Infected Unexposed 
##  34.82208  46.30208

Fecundity

Expectation: infected ones should have many fewer clutches.

# making a plot comparing exposed vs. unexposed
fecundityinfvcontrol <- ggplot(LA_rel_vir_clutches_infvcontrol_long, 
                               aes(x=infvcontrol, y=clutches)) + 
  geom_violin() + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33') +
  xlab("Infected or unexposed") +ylab("Mean # clutches")  +
  theme_cowplot()

fecundityinfvcontrol

LA_rel_vir_clutches_infvcontrol_long$log_clutches <- log((LA_rel_vir_clutches_infvcontrol_long$clutches+1))

infvcontrolrepromodel <- glm(log_clutches ~ infvcontrol, 
             family = "gaussian", 
             data = LA_rel_vir_clutches_infvcontrol_long)

plot(infvcontrolrepromodel) 

summary(infvcontrolrepromodel)
## 
## Call:
## glm(formula = log_clutches ~ infvcontrol, family = "gaussian", 
##     data = LA_rel_vir_clutches_infvcontrol_long)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8405  -0.2869  -0.1015   0.2636   1.5457  
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           0.46921    0.08332   5.632 6.29e-07 ***
## infvcontrolUnexposed  2.06442    0.15726  13.127  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.2846221)
## 
##     Null deviance: 64.702  on 56  degrees of freedom
## Residual deviance: 15.654  on 55  degrees of freedom
## AIC: 94.097
## 
## Number of Fisher Scoring iterations: 2
tapply(LA_rel_vir_clutches_infvcontrol_long$clutches, LA_rel_vir_clutches_infvcontrol_long$infvcontrol, mean)
##   Infected  Unexposed 
##  0.9101239 12.7895833

Let’s arrange the lifespan & fecundity plots into a single figure:

infvcontrolplot <- plot_grid(lifespaninfvcontrol, fecundityinfvcontrol, labels = "auto", ncol = 1, align = "v")
infvcontrolplot

ggsave(here("figures", "infvcontrolplot.pdf"), infvcontrolplot, units = "in", width = 5, height = 8, dpi = 300)
ggsave(here("figures", "infvcontrolplot.jpg"), infvcontrolplot, units = "in", width = 5, height = 8, dpi = 300)

Same time comparison

Let’s first look at the lifespan of animals that were exposed to spores from the same time point

LA_means_sametime <- subset(LA_means, host_time_point == para_time_point)

Lifespan

# plotting animals exposed to parasite from the same time point
lifespansametime <- ggplot(LA_means_sametime, 
                     aes(x=host_time_point, y=mn_lifespan, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  xlab("Host time point") +ylab("Mean lifespan \n(days)") +
  theme_cowplot()

#colors from https://github.com/ciannabp/inauguration/blob/main/R/colors.R 

lifespansametime

mod1a <- glm(mn_lifespan ~ host_time_point, 
             family = "quasipoisson", 
             data = LA_means_sametime)

plot(mod1a) 

summary(mod1a)
## 
## Call:
## glm(formula = mn_lifespan ~ host_time_point, family = "quasipoisson", 
##     data = LA_means_sametime)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6730  -0.3918  -0.1085   0.3645   1.5834  
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.561e+00  7.029e-02  50.661   <2e-16 ***
## host_time_point2  3.363e-02  8.278e-02   0.406    0.688    
## host_time_point3 -7.978e-05  8.038e-02  -0.001    0.999    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.6956699)
## 
##     Null deviance: 16.829  on 26  degrees of freedom
## Residual deviance: 16.573  on 24  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

Lifespan

# plotting animals exposed to parasite from the same time point
rellifespansametime <- ggplot(LA_means_sametime, 
                     aes(x=host_time_point, y=rel_lifespan_mn, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  # ggtitle("Lifespan When Exposed to Parasites from the \nSame Time Point") + 
  xlab("Host time point") +ylab("Relative host lifespan") +
  # labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()

#colors from https://github.com/ciannabp/inauguration/blob/main/R/colors.R 

rellifespansametime

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
mod1arel <- glm(rel_lifespan_mn ~ host_time_point, 
             family = "quasipoisson", 
             data = LA_means_sametime)

plot(mod1arel) 

summary(mod1arel)
## 
## Call:
## glm(formula = rel_lifespan_mn ~ host_time_point, family = "quasipoisson", 
##     data = LA_means_sametime)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.56217  -0.08694  -0.01690   0.10842   0.61409  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)      -0.04584    0.11057  -0.415   0.6821  
## host_time_point2 -0.36181    0.13874  -2.608   0.0154 *
## host_time_point3 -0.17559    0.12927  -1.358   0.1870  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.0467133)
## 
##     Null deviance: 1.4387  on 26  degrees of freedom
## Residual deviance: 1.1079  on 24  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

Reproduction

# plotting animals exposed to parasite from the same time point
reprosametime <- ggplot(LA_means_sametime, 
                     aes(x=host_time_point, y=mn_clutches, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  # ggtitle("Lifespan When Exposed to Parasites from the \nSame Time Point") + 
  xlab("Host time point") +ylab("Mean # clutches") +
  # labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()

reprosametime

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
LA_means_sametime$log_mn_clutches <- log((LA_means_sametime$mn_clutches+1))

mod1b <- glm(log_mn_clutches ~ host_time_point, 
             family = "gaussian", 
             data = LA_means_sametime)

plot(mod1b) 

summary(mod1b)
## 
## Call:
## glm(formula = log_mn_clutches ~ host_time_point, family = "gaussian", 
##     data = LA_means_sametime)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6175  -0.4707  -0.1030   0.1461   1.3115  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)       0.61750    0.31082   1.987   0.0585 .
## host_time_point2  0.08587    0.36776   0.233   0.8174  
## host_time_point3 -0.14681    0.35543  -0.413   0.6832  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.3864272)
## 
##     Null deviance: 9.5874  on 26  degrees of freedom
## Residual deviance: 9.2743  on 24  degrees of freedom
## AIC: 55.771
## 
## Number of Fisher Scoring iterations: 2

Relative reproduction

# plotting animals exposed to parasite from the same time point
relreprosametime <- ggplot(LA_means_sametime, 
                     aes(x=host_time_point, y=rel_clutches_mn, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  xlab("Host time point") +ylab("Relative # clutches") +
  # labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()

relreprosametime

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
LA_means_sametime$log_rel_clutches_mn <- log((LA_means_sametime$rel_clutches_mn+1))

mod1brel <- glm(log_rel_clutches_mn ~ host_time_point, 
             family = "gaussian", 
             data = LA_means_sametime)

plot(mod1brel) 

summary(mod1brel)
## 
## Call:
## glm(formula = log_rel_clutches_mn ~ host_time_point, family = "gaussian", 
##     data = LA_means_sametime)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.12727  -0.07330  -0.04258   0.01503   0.33810  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)       0.12727    0.05986   2.126    0.044 *
## host_time_point2 -0.03734    0.07083  -0.527    0.603  
## host_time_point3 -0.06881    0.06846  -1.005    0.325  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.01433495)
## 
##     Null deviance: 0.35999  on 26  degrees of freedom
## Residual deviance: 0.34404  on 24  degrees of freedom
## AIC: -33.174
## 
## Number of Fisher Scoring iterations: 2

Spores

# plotting animals exposed to parasite from the same time point
sporessametime <- ggplot(LA_means_sametime, 
                     aes(x=host_time_point, y=mn_spores, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  # ggtitle("Lifespan When Exposed to Parasites from the \nSame Time Point") + 
  xlab("Host time point") +ylab("Mean # spores \nper infected host") +
  # labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  scale_y_log10() +
  theme_cowplot()

sporessametime

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
LA_means_sametime$log_mn_spores <- log(LA_means_sametime$mn_spores)

mod1c <- glm(log_mn_spores ~ host_time_point, 
             family = "gaussian", 
             data = LA_means_sametime)

plot(mod1c) 

summary(mod1c)
## 
## Call:
## glm(formula = log_mn_spores ~ host_time_point, family = "gaussian", 
##     data = LA_means_sametime)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.4087  -0.5554   0.3632   0.6789   1.4764  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      11.47023    0.60563  18.939 6.15e-16 ***
## host_time_point2 -0.15376    0.71660  -0.215    0.832    
## host_time_point3 -0.07792    0.69257  -0.113    0.911    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1.467173)
## 
##     Null deviance: 35.287  on 26  degrees of freedom
## Residual deviance: 35.212  on 24  degrees of freedom
## AIC: 91.793
## 
## Number of Fisher Scoring iterations: 2

Parasite growth rate

LA_means_sametime$paragrowth <- (LA_means_sametime$mn_spores)/(LA_means_sametime$mn_lifespan)

paragrowthsametime <- ggplot(LA_means_sametime, 
                     aes(x=host_time_point, y=paragrowth, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  xlab("Host time point") +ylab("Parasite growth rate \n(spores per day)") +
  # labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


paragrowthsametime

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
mod1d <- glm(paragrowth ~ host_time_point, 
             family = "gaussian", 
             data = LA_means_sametime)

plot(mod1d) 

summary(mod1d)
## 
## Call:
## glm(formula = paragrowth ~ host_time_point, family = "gaussian", 
##     data = LA_means_sametime)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4242.4  -1856.4   -411.4   1778.1   6216.9  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)        3719.2     1268.5   2.932  0.00729 **
## host_time_point2    553.5     1500.9   0.369  0.71552   
## host_time_point3   -673.7     1450.6  -0.464  0.64655   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 6436483)
## 
##     Null deviance: 163054339  on 26  degrees of freedom
## Residual deviance: 154475581  on 24  degrees of freedom
## AIC: 504.73
## 
## Number of Fisher Scoring iterations: 2

Combined figure

sametimeplot <- plot_grid(lifespansametime, rellifespansametime, reprosametime, relreprosametime, sporessametime, paragrowthsametime, labels = "auto", ncol = 2, align = "v")
sametimeplot

ggsave(here("figures", "sametimeplot.pdf"), sametimeplot, units = "in", width = 12, height = 7, dpi = 300)
ggsave(here("figures", "sametimeplot.jpg"), sametimeplot, units = "in", width = 12, height = 7, dpi = 300)

Comparing lifespan vs. spore yield

lifespansporeyieldplot <- ggplot(LA_means_sametime, 
                     aes(x=mn_lifespan, y=mn_spores, fill=para_time_point)) + 
  geom_jitter(shape=16, position=position_jitter(width=0.1, height=0), alpha = 0.8, aes(colour=para_time_point)) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  #ggtitle("Spore Yield When Exposed to Parasites from the \nSame Time Point") + 
  xlab("Mean lifespan (days)") +ylab("Mean spore yield \n(spores per host)") +
  labs(color= "Parasite Time Point") + labs(fill="Parasite Time Point") +
  scale_y_log10() +
  geom_smooth(method='lm') +
  theme_cowplot()

lifespansporeyieldplot
## `geom_smooth()` using formula 'y ~ x'

ggsave(here("figures", "lifespansporeyieldplot.pdf"), lifespansporeyieldplot, units = "in", width = 6, height = 3, dpi = 300)
## `geom_smooth()` using formula 'y ~ x'
ggsave(here("figures", "lifespansporeyieldplot.jpg"), lifespansporeyieldplot, units = "in", width = 6, height = 3, dpi = 300)
## `geom_smooth()` using formula 'y ~ x'

Considering tolerance

toleranceplotA <- ggplot(LA_means_sametime, 
                     aes(x=mn_spores, y=mn_lifespan, fill=para_time_point)) + 
  geom_jitter(shape=16, position=position_jitter(width=0.1, height=0), alpha = 0.8, aes(colour=para_time_point)) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  #ggtitle("Spore Yield When Exposed to Parasites from the \nSame Time Point") + 
  xlab("Spore yield") +ylab("Lifespan") +
  labs(color= "Parasite Time Point") + labs(fill="Parasite Time Point") +
  scale_x_log10() +
  geom_smooth(method='lm') +
  theme_cowplot()

#toleranceplotA

#ggsave(here("figures", "lifespansporeyieldplot.pdf"), lifespansporeyieldplot, units = "in", width = 6, height = 3, dpi = 300)
#ggsave(here("figures", "lifespansporeyieldplot.jpg"), lifespansporeyieldplot, units = "in", width = 6, height = 3, dpi = 300)

toleranceplotB <- ggplot(LA_means_sametime, 
                     aes(x=mn_spores, y=mn_clutches, fill=para_time_point)) + 
  geom_jitter(shape=16, position=position_jitter(width=0.1, height=0), alpha = 0.8, aes(colour=para_time_point)) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  #ggtitle("Spore Yield When Exposed to Parasites from the \nSame Time Point") + 
  xlab("Spore yield") +ylab("Clutches") +
  labs(color= "Parasite Time Point") + labs(fill="Parasite Time Point") +
  scale_x_log10() +
  geom_smooth(method='lm') +
  theme_cowplot()

#toleranceplotB

toleranceplot <- plot_grid(toleranceplotA, toleranceplotB, labels = "auto", ncol = 1, align = "v")
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
toleranceplot

LA_means_sametime_time3 <- subset(LA_means_sametime, para_time_point == '3')
LA_means_sametime_time3$log_spores <- log(LA_means_sametime_time3$mn_spores)

cor.test(LA_means_sametime_time3$log_spores, LA_means_sametime_time3$mn_clutches, method="pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  LA_means_sametime_time3$log_spores and LA_means_sametime_time3$mn_clutches
## t = -0.71063, df = 11, p-value = 0.4921
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6817900  0.3860405
## sample estimates:
##        cor 
## -0.2095076
cor.test(LA_means_sametime_time3$log_spores, LA_means_sametime_time3$mn_lifespan, method="pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  LA_means_sametime_time3$log_spores and LA_means_sametime_time3$mn_lifespan
## t = 1.6013, df = 11, p-value = 0.1376
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1527975  0.7952616
## sample estimates:
##       cor 
## 0.4347927

Parasite evolution comparison

Comparing hosts from time 1 and hosts from time 3 that were exposed to parasites from time 3

LA_means_paraevol <- subset(LA_means, host_time_point == 3)

Lifespan

lifespanparaevol <- ggplot(LA_means_paraevol, 
                     aes(x=para_time_point, y=mn_lifespan, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Mean lifespan \n(days)") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


lifespanparaevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
mod2a <- glm(mn_lifespan ~ para_time_point, 
             family = "quasipoisson", 
             data = LA_means_paraevol)

plot(mod2a) 

summary(mod2a)
## 
## Call:
## glm(formula = mn_lifespan ~ para_time_point, family = "quasipoisson", 
##     data = LA_means_paraevol)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.67302  -0.40054  -0.03327   0.30336   1.58336  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.50363    0.03222 108.725   <2e-16 ***
## para_time_point3  0.05734    0.04576   1.253    0.222    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.4831798)
## 
##     Null deviance: 12.774  on 26  degrees of freedom
## Residual deviance: 12.016  on 25  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

Relative lifespan

rellifespanparaevol <- ggplot(LA_means_paraevol, 
                     aes(x=para_time_point, y=rel_lifespan_mn, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Relative host lifespan") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


rellifespanparaevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
mod2arel <- glm(rel_lifespan_mn ~ para_time_point, 
             family = "quasipoisson", 
             data = LA_means_paraevol)

plot(mod2arel) 

summary(mod2arel)
## 
## Call:
## glm(formula = rel_lifespan_mn ~ para_time_point, family = "quasipoisson", 
##     data = LA_means_paraevol)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.25553  -0.05754  -0.01358   0.04401   0.21570  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.28111    0.03458  -8.130 1.75e-08 ***
## para_time_point3  0.05969    0.04908   1.216    0.235    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.01263767)
## 
##     Null deviance: 0.33503  on 26  degrees of freedom
## Residual deviance: 0.31635  on 25  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

Reproduction

reproparaevol <- ggplot(LA_means_paraevol, 
                     aes(x=para_time_point, y=mn_clutches, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Mean # clutches") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


reproparaevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
LA_means_paraevol$log_mn_clutches <- log((LA_means_paraevol$mn_clutches+1))

mod2b <- glm(log_mn_clutches ~ para_time_point, 
             family = "gaussian", 
             data = LA_means_paraevol)

plot(mod2b) 

summary(mod2b)
## 
## Call:
## glm(formula = log_mn_clutches ~ para_time_point, family = "gaussian", 
##     data = LA_means_paraevol)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.47069  -0.25822  -0.03508   0.12837   1.03339  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        0.2582     0.1019   2.535   0.0179 *
## para_time_point3   0.2125     0.1468   1.447   0.1603  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1452936)
## 
##     Null deviance: 3.9366  on 26  degrees of freedom
## Residual deviance: 3.6323  on 25  degrees of freedom
## AIC: 28.462
## 
## Number of Fisher Scoring iterations: 2

Relative reproduction

relreproparaevol <- ggplot(LA_means_paraevol, 
                     aes(x=para_time_point, y=rel_clutches_mn, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Relative # clutches") +
  # labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()

relreproparaevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
LA_means_paraevol$log_rel_clutches_mn <- log((LA_means_paraevol$rel_clutches_mn+1))

mod2brel <- glm(log_rel_clutches_mn ~ para_time_point, 
             family = "gaussian", 
             data = LA_means_paraevol)

plot(mod2brel) 

summary(mod2brel)
## 
## Call:
## glm(formula = log_rel_clutches_mn ~ para_time_point, family = "gaussian", 
##     data = LA_means_paraevol)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.05845  -0.02668  -0.01157   0.01204   0.14514  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)       0.02668    0.01361    1.96   0.0612 .
## para_time_point3  0.03178    0.01961    1.62   0.1178  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.002593259)
## 
##     Null deviance: 0.071638  on 26  degrees of freedom
## Residual deviance: 0.064831  on 25  degrees of freedom
## AIC: -80.236
## 
## Number of Fisher Scoring iterations: 2

Spores

sporesparaevol <- ggplot(LA_means_paraevol, 
                     aes(x=para_time_point, y=mn_spores, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Mean # spores \nper infected host") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  scale_y_log10() +
  theme_cowplot()


sporesparaevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
LA_means_paraevol$log_mn_spores <- log(LA_means_paraevol$mn_spores)

mod2c <- glm(log_mn_spores ~ para_time_point, 
             family = "gaussian", 
             data = LA_means_paraevol)

plot(mod2c) 

summary(mod2c)
## 
## Call:
## glm(formula = log_mn_spores ~ para_time_point, family = "gaussian", 
##     data = LA_means_paraevol)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1748  -0.4538   0.1635   0.4625   0.8994  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       11.8023     0.1643  71.817   <2e-16 ***
## para_time_point3  -0.4099     0.2368  -1.731   0.0958 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.3780972)
## 
##     Null deviance: 10.5853  on 26  degrees of freedom
## Residual deviance:  9.4524  on 25  degrees of freedom
## AIC: 54.284
## 
## Number of Fisher Scoring iterations: 2

Parasite growth rate

LA_means_paraevol$paragrowth <- (LA_means_paraevol$mn_spores)/(LA_means_paraevol$mn_lifespan)

paragrowthparaevol <- ggplot(LA_means_paraevol, 
                     aes(x=para_time_point, y=paragrowth, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Parasite growth rate \n(spores per day)") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


paragrowthparaevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
mod2d <- glm(paragrowth ~ para_time_point, 
             family = "gaussian", 
             data = LA_means_paraevol)

plot(mod2d) 

summary(mod2d)
## 
## Call:
## glm(formula = paragrowth ~ para_time_point, family = "gaussian", 
##     data = LA_means_paraevol)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3060.5  -1444.9    126.3   1405.7   4602.6  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        4482.9      511.7   8.762  4.3e-09 ***
## para_time_point3  -1437.4      737.4  -1.949   0.0626 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 3665074)
## 
##     Null deviance: 105553897  on 26  degrees of freedom
## Residual deviance:  91626856  on 25  degrees of freedom
## AIC: 488.63
## 
## Number of Fisher Scoring iterations: 2

Combined figure

paraevolplot <- plot_grid(lifespanparaevol, rellifespanparaevol, reproparaevol, relreproparaevol, sporesparaevol, paragrowthparaevol, labels = "auto", ncol = 2, align = "v")
paraevolplot

ggsave(here("figures", "paraevolplot.pdf"), paraevolplot, units = "in", width = 10, height = 7, dpi = 300)
ggsave(here("figures", "paraevolplot.jpg"), paraevolplot, units = "in", width = 10, height = 7, dpi = 300)

Host evolution comparison

Comparing hosts from time 1 and hosts from time 3 that were exposed to parasites from time 3

LA_means_hostevol <- subset(LA_means, para_time_point == 1)

Lifespan

lifespanhostevol <- ggplot(LA_means_hostevol, 
                     aes(x=host_time_point, y=mn_lifespan, fill=host_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#808080")) + 
  scale_color_manual(values=c("#5445b1", "#808080")) +
  xlab("Host time point") +ylab("Mean lifespan \n(days)") +
  labs(color= "Host time point") + labs(fill="Host time point") +
  theme_cowplot()


lifespanhostevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
mod3a <- glm(mn_lifespan ~ host_time_point, 
             family = "quasipoisson", 
             data = LA_means_hostevol)

plot(mod3a) 

summary(mod3a)
## 
## Call:
## glm(formula = mn_lifespan ~ host_time_point, family = "quasipoisson", 
##     data = LA_means_hostevol)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.41726  -0.38792   0.03013   0.35122   1.11192  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.56105    0.05437  65.497   <2e-16 ***
## host_time_point3 -0.05742    0.06205  -0.925    0.369    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.4162173)
## 
##     Null deviance: 7.1131  on 17  degrees of freedom
## Residual deviance: 6.7604  on 16  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

Relative lifespan

rellifespanhostevol <- ggplot(LA_means_hostevol, 
                     aes(x=host_time_point, y=rel_lifespan_mn, fill=host_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#808080")) + 
  scale_color_manual(values=c("#5445b1", "#808080")) +
  xlab("Host time point") +ylab("Relative host lifespan") +
  labs(color= "Host time point") + labs(fill="Host time point") +
  theme_cowplot()


rellifespanhostevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
mod3arel <- glm(rel_lifespan_mn ~ host_time_point, 
             family = "quasipoisson", 
             data = LA_means_hostevol)

plot(mod3arel) 

summary(mod3arel)
## 
## Call:
## glm(formula = rel_lifespan_mn ~ host_time_point, family = "quasipoisson", 
##     data = LA_means_hostevol)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.56217  -0.05907  -0.01070   0.04401   0.61409  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## (Intercept)      -0.04584    0.11770  -0.389    0.702
## host_time_point3 -0.23527    0.13733  -1.713    0.106
## 
## (Dispersion parameter for quasipoisson family taken to be 0.05292627)
## 
##     Null deviance: 0.98270  on 17  degrees of freedom
## Residual deviance: 0.83318  on 16  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

Reproduction

reprohostevol <- ggplot(LA_means_hostevol, 
                     aes(x=host_time_point, y=mn_clutches, fill=host_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#808080")) + 
  scale_color_manual(values=c("#5445b1", "#808080")) +
  xlab("Host time point") +ylab("Mean # clutches") +
  labs(color= "Host time point") + labs(fill="Host time point") +
  theme_cowplot()


reprohostevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
LA_means_hostevol$log_mn_clutches <- log((LA_means_hostevol$mn_clutches+1))

mod3b <- glm(log_mn_clutches ~ host_time_point, 
             family = "gaussian", 
             data = LA_means_hostevol)

plot(mod3b) 

summary(mod3b)
## 
## Call:
## glm(formula = log_mn_clutches ~ host_time_point, family = "gaussian", 
##     data = LA_means_hostevol)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.61750  -0.25822  -0.03508   0.10674   1.22833  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        0.6175     0.2192   2.817   0.0124 *
## host_time_point3  -0.3593     0.2485  -1.446   0.1676  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1921489)
## 
##     Null deviance: 3.4760  on 17  degrees of freedom
## Residual deviance: 3.0744  on 16  degrees of freedom
## AIC: 25.271
## 
## Number of Fisher Scoring iterations: 2

Relative reproduction

relreprohostevol <- ggplot(LA_means_hostevol, 
                     aes(x=host_time_point, y=rel_clutches_mn, fill=host_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#808080")) + 
  scale_color_manual(values=c("#5445b1", "#808080")) +
  xlab("Host time point") +ylab("Relative # clutches") +
  labs(color= "Host time point") + labs(fill="Host time point") +
  theme_cowplot()


relreprohostevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
LA_means_hostevol$log_rel_clutches_mn <- log((LA_means_hostevol$rel_clutches_mn+1))

mod3brel <- glm(log_rel_clutches_mn ~ host_time_point, 
             family = "gaussian", 
             data = LA_means_hostevol)

plot(mod3brel) 

summary(mod3brel)
## 
## Call:
## glm(formula = log_rel_clutches_mn ~ host_time_point, family = "gaussian", 
##     data = LA_means_hostevol)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.12727  -0.02668  -0.00902   0.01072   0.33810  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)       0.12727    0.05063   2.513   0.0230 *
## host_time_point3 -0.10059    0.05741  -1.752   0.0989 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.01025551)
## 
##     Null deviance: 0.19557  on 17  degrees of freedom
## Residual deviance: 0.16409  on 16  degrees of freedom
## AIC: -27.477
## 
## Number of Fisher Scoring iterations: 2

Spores

sporeshostevol <- ggplot(LA_means_hostevol, 
                     aes(x=host_time_point, y=mn_spores, fill=host_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#808080")) + 
  scale_color_manual(values=c("#5445b1", "#808080")) +
  xlab("Host time point") +ylab("Mean # spores \nper infected host") +
  labs(color= "Host time point") + labs(fill="Host time point") +
  scale_y_log10() +
  theme_cowplot()


sporeshostevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
LA_means_hostevol$log_mn_spores <- log(LA_means_hostevol$mn_spores)

mod3c <- glm(log_mn_spores ~ host_time_point, 
             family = "gaussian", 
             data = LA_means_hostevol)

plot(mod3c) 

summary(mod3c)
## 
## Call:
## glm(formula = log_mn_spores ~ host_time_point, family = "gaussian", 
##     data = LA_means_hostevol)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1748  -0.3479   0.1552   0.4171   1.0125  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        11.470      0.306  37.480   <2e-16 ***
## host_time_point3    0.332      0.347   0.957    0.353    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.3746263)
## 
##     Null deviance: 6.337  on 17  degrees of freedom
## Residual deviance: 5.994  on 16  degrees of freedom
## AIC: 37.289
## 
## Number of Fisher Scoring iterations: 2

Parasite growth rate

LA_means_hostevol$paragrowth <- (LA_means_hostevol$mn_spores)/(LA_means_hostevol$mn_lifespan)

paragrowthhostevol <- ggplot(LA_means_hostevol, 
                     aes(x=host_time_point, y=paragrowth, fill=host_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#808080")) + 
  scale_color_manual(values=c("#5445b1", "#808080")) +
  xlab("Host time point") +ylab("Parasite growth rate \n(spores per day)") +
  labs(color= "Host time point") + labs(fill="Host time point") +
  theme_cowplot()


paragrowthhostevol

# ggsave(here("figures", "lifespansametime.pdf"), lifespansametime, units = "in", width = 7, height = 5, dpi = 300)
mod3d <- glm(paragrowth ~ host_time_point, 
             family = "gaussian", 
             data = LA_means_hostevol)

plot(mod3d) 

summary(mod3d)
## 
## Call:
## glm(formula = paragrowth ~ host_time_point, family = "gaussian", 
##     data = LA_means_hostevol)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3060.5  -1781.2    163.5   1234.8   4602.6  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)        3719.2     1119.4   3.322  0.00431 **
## host_time_point3    763.7     1269.3   0.602  0.55581   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 5012430)
## 
##     Null deviance: 82013591  on 17  degrees of freedom
## Residual deviance: 80198882  on 16  degrees of freedom
## AIC: 332.66
## 
## Number of Fisher Scoring iterations: 2

Combined figure

hostevolplot <- plot_grid(lifespanhostevol, rellifespanhostevol, reprohostevol, relreprohostevol, sporeshostevol, paragrowthhostevol, labels = "auto", ncol = 2, align = "v")
hostevolplot

ggsave(here("figures", "hostevolplot.pdf"), hostevolplot, units = "in", width = 10, height = 7, dpi = 300)
ggsave(here("figures", "hostevolplot.jpg"), hostevolplot, units = "in", width = 10, height = 7, dpi = 300)